% Figure 2 in "Commodity Prices, Convenience Yields and Inflation"
% March 2011

clear all;

load cyj.dat;                   % convenience yields on 23 commodities
load us_intrate_monthly.dat;    % US 3-month T-bill rate
load imf_index_monthly.dat;     % IMF aggregate commodity price index
load usd_monthly_ave.dat;       % USD trade-weighted index

i=us_intrate_monthly;
SA=imf_index_monthly;
er=usd_monthly_ave;
nn=rows(er);

h=1;            % forecast horizon
rhat=2;         % number of principal components
nwl=4;          % number of lags in Newey-West HAC estimator

[ehat,Fhat,lamhat,ve]=pc(standard(cyj),rhat); % PC estimation

ds_h=100*(log(SA(3+h:nn,1))-log(SA(3:nn-h,1)));
ds_1=100*(log(SA(3:nn-h,1))-log(SA(2:nn-h-1,1)));
ds_2=100*(log(SA(2:nn-h-1,1))-log(SA(1:nn-h-2,1)));
dx=100*(log(er(3:nn-h,1))-log(er(2:nn-h-1,1)));
ir=i(3:nn-h,1);
p1_cy=Fhat(3:nn-h,1);
p2_cy=Fhat(3:nn-h,2);
nobs=rows(ds_h);

% Commodity price equation
x=[p1_cy p2_cy ir dx ds_1 ds_2];
res=nwest(ds_h,[ones(nobs,1) x],nwl);
x1=[ones(19,1) p1_cy(nobs-18)*ones(19,1) p2_cy(nobs-18)*ones(19,1) ir(nobs-18:nobs) dx(nobs-18:nobs) ds_1(nobs-18:nobs) ds_2(nobs-18:nobs)];
x2=[ones(19,1) p1_cy(nobs-18:nobs) p2_cy(nobs-18:nobs) ir(nobs-18)*ones(19,1) dx(nobs-18:nobs) ds_1(nobs-18:nobs) ds_2(nobs-18:nobs)];
x3=[ones(19,1) p1_cy(nobs-18:nobs) p2_cy(nobs-18:nobs) ir(nobs-18:nobs) dx(nobs-18)*ones(19,1) ds_1(nobs-18:nobs) ds_2(nobs-18:nobs)];
x4=[ones(19,1) p1_cy(nobs-18)*ones(19,1) p2_cy(nobs-18)*ones(19,1) ir(nobs-18)*ones(19,1) dx(nobs-18)*ones(19,1) ds_1(nobs-18:nobs) ds_2(nobs-18:nobs)];
yhat1=x1*res.beta;
yhat2=x2*res.beta;
yhat3=x3*res.beta;
yhat4=x4*res.beta;

act=[]; cf1=[]; cf2=[]; cf3=[]; cf4=[];
act(1,:)=SA(rows(SA)-18,:);
cf1(1,:)=SA(rows(SA)-18,:);
cf2(1,:)=SA(rows(SA)-18,:);
cf3(1,:)=SA(rows(SA)-18,:);
cf4(1,:)=SA(rows(SA)-18,:);
for j=2:20;
    act(j,:)=act(j-1,:)*(1+ds_h(nobs-20+j,:)/100);
    cf1(j,:)=cf1(j-1,:)*(1+yhat1(j-1,:)/100);
    cf2(j,:)=cf2(j-1,:)*(1+yhat2(j-1,:)/100);
    cf3(j,:)=cf3(j-1,:)*(1+yhat3(j-1,:)/100);
    cf4(j,:)=cf4(j-1,:)*(1+yhat4(j-1,:)/100);
end;

xaxis=seqa(1,1,rows(cf1));
plot(act(2:20,1),'-k','LineWidth',1)
hold on
plot(cf1(2:20,1),'-bs','MarkerSize',4);
plot(cf2(2:20,1),'-go','MarkerSize',4);
plot(cf3(2:20,1),'-r*','MarkerSize',4);
plot(cf4(2:20,1),'-cd','MarkerSize',4);
hold off
set(gca,'XTick',[1 6 12 18])
set(gca,'XTickLabel',{'Jan 2007' 'June 2007' ...
    'Dec 2007' 'June 2008'})
h = legend('actual values of commodity price index','counterfactual: fixed convenience yield',...,
    'counterfactual: fixed interest rate','counterfactual: fixed exchange rate',...,
    'counterfactual: all fixed',2);
set(h,'Interpreter','none')
